REFMAP Listening test 1 Part A analysis: Random forest variable importance identification (median responses)

This notebook analyses the Part A data in terms of variable importance using a random forest model based on conditional inference trees and a conditional permutation variable importance algorithm.

Setup

Load packages

# load packages
library(tidyr)
library(ggplot2)
library(party)
library(conflicted)
library(tidyverse)
library(openxlsx)
library(caret)
library(viridis)
library(cowplot)
# set package parameters
theme_set(theme_bw())

# plot colour scheme

mycolourlist = list(c(0, 102, 255), c(0, 204, 153), c(255, 0, 102), c(74, 111, 152), c(251, 164, 49), c(204, 153, 255), c(90, 192, 255), c(80, 245, 233), c(255, 90, 192), c(164, 201, 242), c(255, 254, 139), c(255, 243, 255))
mycolours = matrix()

for (ii in 1:length(mycolourlist)){
  mycolours[ii] = rgb(mycolourlist[[ii]][1]/255,
                      mycolourlist[[ii]][2]/255,
                      mycolourlist[[ii]][3]/255)
}

# toggle to save plots
saveplots = TRUE

if (saveplots){
  # set output plot directory
  choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2))  # workaround for bug in RTerm choose.dir
  outFigPath <- utils::choose.dir(caption="Select output folder to save plots '03 Experiment\\Experiment 1\\Analysis\\Plots'")
  
  if (!dir.exists(file.path(outFigPath, "svg"))){dir.create(file.path(outFigPath, "svg"))}
  if (!dir.exists(file.path(outFigPath, "pdf"))){dir.create(file.path(outFigPath, "pdf"))}
  
}

# toggle to save data
savedata = TRUE

if (savedata){
  # set output plot directory
  if (saveplots==FALSE){
    choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2))  # workaround for bug in RTerm choose.dir
  }
  outDataPath <- utils::choose.dir(caption="Select output folder to save data '03 Experiment\\Experiment 1\\Analysis\\R'")
}
 

Import data and wrangle

Part A

stimDataApath <- utils::choose.files(caption=r"(Select refmap_listest1_testdataA_ByStim.csv from 03 Experiment\Experiment 1\Analysis\PostProcess)",
                                     filters=matrix(data=c("refmap_listest1_testdataA_ByStim.csv", "refmap_listest1_testdataA_ByStim.csv"), ncol=2))

stimNoticeDataApath <- utils::choose.files(caption=r"(Select refmap_listest1_testdataANoticeFilt_ByStim.csv from 03 Experiment\Experiment 1\Analysis\PostProcess)", filters=matrix(data=c("refmap_listest1_testdataANoticeFilt_ByStim.csv", "refmap_listest1_testdataANoticeFilt_ByStim.csv"), ncol=2))

stimDataA <- utils::read.csv(stimDataApath, header=TRUE)
stimNoticeDataA <- utils::read.csv(stimNoticeDataApath, header=TRUE)

colnames(stimDataA)[1] <- "Stimulus"
colnames(stimNoticeDataA)[1] <- "Stimulus"

# make response proportions into percentages
stimDataA[['HighAnnoyPc']] <- stimDataA[['HighAnnoyProp']]*100
stimDataA[['dHighAnnoyPc']] <- stimDataA[['dHighAnnoyProp']]*100

# make response proportions into percentages
stimNoticeDataA[['HighAnnoyPcFilt']] <- stimNoticeDataA[['HighAnnoyPropFilt']]*100
stimNoticeDataA[['dHighAnnoyPcFilt']] <- stimNoticeDataA[['dHighAnnoyPropFilt']]*100
stimNoticeDataA[['NoticedPcFilt']] <- stimNoticeDataA[['NoticedPropFilt']]*100

# function to encode categorical to ordinal numeric variables
encode_ordinal <- function(x, order=unique(x)) {
  x <- as.numeric(factor(x, levels=order, exclude=NULL, order=TRUE))
  x
}

# definition of ordinal variable levels
SNRCatsA <- c("No UAS", "-16", "-10", "-4", "2", "8")
UASLAeqCatsA <- c("No UAS", "42", "48", "54", "60")

The aggregated data by stimulus are assigned to a dataframe, relevant categorical variables are converted to ordinal, and then the variable subset of interest is selected, NA rows dropped (ie, the ‘no UAS’ stimuli, as the conditional variable importance algorithm cannot currently handle NA values, which are present in all the UAS dB metrics), and a formula assigned.

stimDataANum <- data.frame()

stimDataANum <- cbind(stimDataA[, 'Stimulus'],
                      stimDataA[, which(colnames(stimDataA)=="UASLAeq"):
                                  which(colnames(stimDataA)=="SNRlevel")],
                      stimDataA[, "IntermitRatioMaxLR", drop=F],
                      stimDataA[, which(colnames(stimDataA)=="UASLAeqMaxLR"):
                                  which(colnames(stimDataA)=="UASImpulsSHM05ExMaxLR")],
                      stimDataA[, which(colnames(stimDataA)=="LAeqLAF90diff"):
                                  which(colnames(stimDataA)=="dImpulsSHM05ExMaxLR")],
                      stimDataA[, which(colnames(stimDataA)=="ValenceMedian"):
                                  which(colnames(stimDataA)=="dHighAnnoyProp")],
                      stimDataA[, which(colnames(stimDataA)=="HighAnnoyPc"):
                                  which(colnames(stimDataA)=="dHighAnnoyPc")])

# remove duplicated variables
stimDataANum <- subset(stimDataANum, select = -c(UASLAeqMaxLR, UASLAEMaxLR))

colnames(stimDataANum)[1] <- "Stimulus"

# make the categorical outcome variables numerical dummy factors
# stimDataANum[['SNRlevel']] <- encode_ordinal(stimDataA[['SNRlevel']], order=SNRCatsA)
# stimDataANum[['UASLAeq']] <- encode_ordinal(stimDataA[['UASLAeq']], order=UASLAeqCatsA)

# make the discrete ordinal outcome variables factors
stimDataANum[['ValenceMedian']] <- factor(stimDataANum[['ValenceMedian']], levels=c(1, 2, 3, 4, 5), order=TRUE)
stimDataANum[['ArousalMedian']] <- factor(stimDataANum[['ArousalMedian']], levels=c(1, 2, 3, 4, 5), order=TRUE)
stimDataANum[['AnnoyMedian']] <- factor(stimDataANum[['AnnoyMedian']], levels=c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), order=TRUE)
stimDataANum[['dValenceMedian']] <- factor(stimDataANum[['dValenceMedian']], levels=c(-4, -3, -2, -1, 0, 1, 2, 3, 4), order=TRUE)
stimDataANum[['dArousalMedian']] <- factor(stimDataANum[['dArousalMedian']], levels=c(-4, -3, -2, -1, 0, 1, 2, 3, 4), order=TRUE)
stimDataANum[['dAnnoyMedian']] <- factor(stimDataANum[['dAnnoyMedian']], levels=c(-10, -9, -8, -7, -6, -5, -4, -3, -2, -1,
                                                                                 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), order=TRUE)

stimNoticeDataANum <- data.frame()

stimNoticeDataANum <- cbind(stimNoticeDataA[, 'Stimulus'],
                            stimNoticeDataA[, which(colnames(stimNoticeDataA)=="UASLAeq"):
                                              which(colnames(stimNoticeDataA)=="SNRlevel")],
                            stimNoticeDataA[, "IntermitRatioMaxLR", drop=F],
                            stimNoticeDataA[, which(colnames(stimNoticeDataA)=="UASLAeqMaxLR"):
                                              which(colnames(stimNoticeDataA)=="UASImpulsSHM05ExMaxLR")],
                            stimNoticeDataA[, which(colnames(stimNoticeDataA)=="LAeqLAF90diff"):
                                              which(colnames(stimNoticeDataA)=="dImpulsSHM05ExMaxLR")],
                            stimNoticeDataA[, which(colnames(stimNoticeDataA)=="NoticedTotalFilt"):
                                              which(colnames(stimNoticeDataA)=="NoticedPcFilt")])

stimNoticeDataANum <- subset(stimNoticeDataANum, select = -c(UASLAeqMaxLR, UASLAEMaxLR))

colnames(stimNoticeDataANum)[1] <- "Stimulus"

# make the categorical outcome variables numerical dummy factors
# stimNoticeDataANum[['SNRlevel']] <- encode_ordinal(stimNoticeDataA[['SNRlevel']], order=SNRCatsA)
# stimNoticeDataANum[['UASLAeq']] <- encode_ordinal(stimNoticeDataA[['UASLAeq']], order=UASLAeqCatsA)

Random forest functions

Write a function to train a conditional-inference random forest (crf) model on input data according to input formula, iterate over input random seeds, average error and variable importance metrics, and output metrics with plotted

Averaging over multiple random seeds

multi_crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
  # initialise variables
  crfOOBErrAll <- 0
  crfOOBRMSE <- 0
  crfOOBMAE <- 0
  crfOOBErrR2 <- 0
  crfMarPermImpVals <- 0
  crfConPermImpVals <- 0
  crfMarPermImpValsPerTree <- data.frame()
  crfConPermImpValsPerTree <- data.frame()
  
  for (iters in 1:length(seeds)){
    
    # formula for regression
    formVars <- reformulate(iVars, dVar)
    
    # set random seed
    set.seed(seeds[iters])
    # train crf model
    crfModel <- party::cforest(formVars, data=dataIn,
                               controls=party::cforest_unbiased(ntree=ntree,
                                                                mtry=mtry,
                                                                minsplit=minsplit,
                                                                minbucket=minbucket))
    
    # get OOB predictions
    crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
    
    # get OOB error
    crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
                                            - as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))

    # OOB RMSE, error quartiles and Rsquared
    crfOOBRMSE <- crfOOBRMSE + sqrt(mean(crfModelOOBErr^2))
    crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
    crfOOBErrR2 <- crfOOBErrR2 + cor(as.numeric(as.matrix(crfModelOOB)),
                                     as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2

    # set random seed
    set.seed(seeds[iters])

    # set random seed
    set.seed(seeds[iters])
    # conditional variable permutation importance
    crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
    
    crfConPermImpVals <- crfConPermImpVals + crfConPermImp$values
    crfConPermImpValsPerTree <- rbind(crfConPermImpValsPerTree, crfConPermImp$perTree)
  }
  
  # average metrics
  crfOOBErrAll <- crfOOBErrAll/length(seeds)
  crfOOBRMSE <- crfOOBRMSE/length(seeds)
  crfOOBMAE <- crfOOBMAE/length(seeds)
  crfOOBErrR2 <- crfOOBErrR2/length(seeds)
  crfConPermImpVals <- data.frame(CondPermImp=sort(crfConPermImpVals/length(seeds), decreasing=TRUE))
  #crfConPermImpValsMn <- sort(colMeans(crfConPermImpValsPerTree), decreasing=TRUE)  # mean of conditional variable importance values per tree (same as crfConPermImpVals)
  crfConPermImpValsQtl <- data.frame(apply(crfConPermImpValsPerTree, 2, quantile, probs=c(0.25, 0.50, 0.75)))
  
  resultsOut <- list('OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals,                      'conditional_permimp_perTree'=crfConPermImpValsPerTree, 'conditional_permimp_qtl'=crfConPermImpValsQtl)
  return(resultsOut)
}

Comparing rankings from two seeds

crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
  # initialise variables
  crfOOBErrAll <- 0
  crfOOBRMSE <- 0
  crfOOBMAE <- 0
  crfOOBErrR2 <- 0
  crfMarPermImpVals <- 0
  crfConPermImpVals <- 0
  crfMarPermImpValsPerTree <- data.frame()
  crfConPermImpValsPerTree <- data.frame()

  # formula for regression
  formVars <- reformulate(iVars, dVar)
  
  for (iters in 1:length(seeds)){
  
    # set random seed
    set.seed(seeds[iters])
    # train crf model
    crfModel <- party::cforest(formVars, data=dataIn,
                               controls=party::cforest_unbiased(ntree=ntree,
                                                                mtry=mtry,
                                                                minsplit=minsplit,
                                                                minbucket=minbucket))
    
    # conditional variable permutation importance
    crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
    
    crfConPermImpVals <- crfConPermImp$values
    
    if (iters == 1){
      crfConPermImpVals1 <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
      crfConPermImpValsPerTree1 <- crfConPermImp$perTree
      crfConPermImpValsQtl1 <- data.frame(apply(crfConPermImpValsPerTree1, 2, quantile, probs=c(0.25, 0.50, 0.75)))
      
      # get OOB predictions
      crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
      
      # get OOB error
      crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
                                              - as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))
      
      # OOB RMSE, error quartiles and Rsquared
      crfOOBRMSE <- sqrt(mean(crfModelOOBErr^2))
      crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
      crfOOBErrR2 <- cor(as.numeric(as.matrix(crfModelOOB)),
                                    as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2

      }
    
    else{
      crfConPermImpValsN <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
      
      nVarImpChecks <- min(max(sum(crfConPermImpVals1 >= crfConPermImpVals1$CondPermImp[1]*0.1),
                               sum(crfConPermImpValsN >= crfConPermImpValsN$CondPermImp[1]*0.1)), 4)  # number of variable importance values with a value at least 10% of the highest importance
      if (sum(rownames(crfConPermImpVals1)[1:nVarImpChecks] != rownames(crfConPermImpValsN)[1:nVarImpChecks]) > 0){
        warning("Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')")
      }
      else{
        resultsOut <- list('OOB_errors'=crfModelOOBErr, 'OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals1, 'conditional_permimp_perTree'=crfConPermImpValsPerTree1, 'conditional_permpimp_qtl'=crfConPermImpValsQtl1)
        return(resultsOut)
      }
      
    }
    
  }

}

Hyperparameter tuning

mtryTune <- function(dataIn, iVars, dVar, seed, ntrees, minsplit=20, minbucket=7){
  
  formVars <- reformulate(iVars, dVar)
  
  # set mtry values and corresponding iVars/mtry ratios
  iVars_mtrys <- c(10.5, 5.25, 3.5, 2.75, 2.25, 1.75, 1.5, 1.25)
  mtrys <- as.integer(length(iVars)/iVars_mtrys)
  iVars_mtrys <- iVars_mtrys[mtrys >= 2]  # remove 0 or 1 values
  mtrys <- mtrys[mtrys >= 2]  # remove 0 or 1 values
  
  # remove any duplicated values
  iVars_mtrys <- iVars_mtrys[!(duplicated(mtrys) | duplicated(mtrys, fromLast = TRUE))]
  mtrys <- mtrys[!(duplicated(mtrys) | duplicated(mtrys, fromLast = TRUE))]

  # ensure mtrys is less than length(iVars)
  iVars_mtrys <- iVars_mtrys[mtrys < length(iVars)]
  mtrys <- mtrys[mtrys < length(iVars)]

  resRMSEMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
  resRsquaredMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
  resMAEMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
  
  for (ii in seq(1, length(ntrees))){
    set.seed(seed)
    ntree = ntrees[ii]
    tuneMod <- caret::train(formVars,
                            data=dataIn,
                            method='cforest',
                            controls=party::cforest_unbiased(ntree=ntree,
                                                             minsplit=minsplit,
                                                             minbucket=minbucket),
                            tuneGrid=data.frame(.mtry=mtrys),
                            trControl = trainControl(method = "oob"))
    
    tuneMod$results <- cbind(tuneMod$results, data.frame(iVars_mtrys=iVars_mtrys))
    
    print(tuneMod$results)
    print(tuneMod$bestTune)
    
    resRMSEMap[, ii] = tuneMod$results$RMSE
    resRsquaredMap[, ii] = tuneMod$results$Rsquared
    resMAEMap[, ii] = tuneMod$results$MAE
    
  }
  
  
  # convert to data frames with mtry as row names and ntree as column names, and convert to long format using tidyverse
  resdfRMSEMap <- as.data.frame(resRMSEMap)
  rownames(resdfRMSEMap) <- mtrys
  colnames(resdfRMSEMap) <- ntrees
  resdfRsquaredMap <- as.data.frame(resRsquaredMap)
  rownames(resdfRsquaredMap) <- mtrys
  colnames(resdfRsquaredMap) <- ntrees
  resdfMAEMap <- as.data.frame(resMAEMap)
  rownames(resdfMAEMap) <- mtrys
  colnames(resdfMAEMap) <- ntrees
  
  
  # convert dataframes to long format using tidyverse
  resdfRMSEMap <- resdfRMSEMap |>
                      rownames_to_column('mtry') |>
                          gather(key='ntree', value='RMSE', -mtry)
  
  resdfRsquaredMap <- resdfRsquaredMap |>
                          rownames_to_column('mtry') |>
                              gather(key='ntree', value='Rsquared', -mtry)
  
  resdfMAEMap <- resdfMAEMap |>
                    rownames_to_column('mtry') |>
                        gather(key='ntree', value='MAE', -mtry)
  
  # ensure ntree and mtry columns are ordered factors
  resdfRMSEMap$ntree <- factor(resdfRMSEMap$ntree, levels=as.character(ntrees))
  resdfRMSEMap$mtry <- factor(resdfRMSEMap$mtry, levels=as.character(mtrys))
  
  resdfRsquaredMap$ntree <- factor(resdfRsquaredMap$ntree, levels=as.character(ntrees))
  resdfRsquaredMap$mtry <- factor(resdfRsquaredMap$mtry, levels=as.character(mtrys))
  
  resdfMAEMap$ntree <- factor(resdfMAEMap$ntree, levels=as.character(ntrees))
  resdfMAEMap$mtry <- factor(resdfMAEMap$mtry, levels=as.character(mtrys))
  
  # plot heatmaps using ggplot, with extreme (min or max) value plotted as overlaid point using annotate and colourbar scale reversed
  pHeatmapRMSE <- ggplot(resdfRMSEMap) +
                    geom_tile(aes(x=ntree, y=mtry, fill=RMSE)) +
                        scale_fill_viridis(option="viridis", direction=-1) +
                          geom_point(data=resdfRMSEMap[which(resdfRMSEMap$RMSE
                                                             == min(resdfRMSEMap$RMSE),
                                                             arr.ind = TRUE),],
                                     aes(x=ntree, y=mtry), colour="red", size=2) +
                            guides(colour = guide_colourbar(reverse=TRUE)) +
                              labs(x="ntree", y="mtry", fill="RMSE") +
                                theme(text = element_text(family = "serif"))
  
  pHeatmapRsquared <- ggplot(resdfRsquaredMap) +
                        geom_tile(aes(x=ntree, y=mtry, fill=Rsquared)) +
                            scale_fill_viridis(option="viridis", direction=1) +
                              geom_point(data=resdfRsquaredMap[which(resdfRsquaredMap$Rsquared
                                                                     == max(resdfRsquaredMap$Rsquared),
                                                                     arr.ind = TRUE),],
                                         aes(x=ntree, y=mtry), colour="red", size=2) +
                                guides(colour = guide_colourbar(reverse=TRUE)) +
                                  labs(x="ntree", y="mtry", fill="Rsquared") +
                                    theme(text = element_text(family = "serif"))
  
  pHeatmapMAE <- ggplot(resdfMAEMap) +
                    geom_tile(aes(x=ntree, y=mtry, fill=MAE)) +
                        scale_fill_viridis(option="viridis", direction=-1) +
                          geom_point(data=resdfMAEMap[which(resdfMAEMap$MAE
                                                            == min(resdfMAEMap$MAE),
                                                            arr.ind = TRUE),],
                                     aes(x=ntree, y=mtry), colour="red", size=2) +
                            guides(colour = guide_colourbar(reverse=TRUE)) +
                              labs(x="ntree", y="mtry", fill="MAE") +
                                theme(text = element_text(family = "serif"))
  
  p <-  cowplot::plot_grid(pHeatmapRMSE,
                           pHeatmapRsquared,
                           pHeatmapMAE,
                           ncol=3, nrow=1)
  
  return(p)
  
}  # end of function

Part A analysis

Set global parameters

permImpCondThres <- 0.95
minsplit <- 20
minbucket <- 7
ntrees <- c(251, 501, 1001, 1501, 2501, 4001, 5501)

Mean change in annoyance

Initialise results output variables

resdAnnoyMnFit <- data.frame(RMSE = numeric(),
                             MAE = numeric(),
                             Rsquared = numeric())
resdAnnoyMnPermImp <- list()

All variables (absolute and difference)

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)

Set variables

iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
           names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dAnnoyMean"

seeds <- c(14569, 98651, 54654498, 454948, 41321)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAllVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAllVarsHyperTune.svg")

  ggsave(filename="PtAdAnnoyMnAllVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAllVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 10

resultsOutAbsDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 0.434307

resultsOutAbsDiffs$OOB_MAE

[1] 0.3533671

resultsOutAbsDiffs$Rsquared

[1] 0.8738479

Train multiple seeds model

resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 0.430126

resultsOutAbsDiffs$OOB_MAE

[1] 0.3537057

resultsOutAbsDiffs$Rsquared

[1] 0.875829

# store results
resdAnnoyMnFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdAnnoyMnFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdAnnoyMnFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdAnnoyMnPermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAllVarsConPermimp.svg")

  ggsave(filename="PtAdAnnoyMnAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAllVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAllVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAllVarsConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnAllVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAllVarsConPermimp.pdf")
}

Absolute variables

Set variables

iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "dAnnoyMean"

seeds <- c(578312, 544, 84894, 54654, 153157)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsVarsHyperTune.svg")

  ggsave(filename="PtAdAnnoyMnAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutAbs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 0.5920291

resultsOutAbs$OOB_MAE

[1] 0.4797379

resultsOutAbs$Rsquared

[1] 0.8046045

Train multiple seeds model

resultsOutAbs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 0.5898948

resultsOutAbs$OOB_MAE

[1] 0.4771605

resultsOutAbs$Rsquared

[1] 0.8064115

# store results
resdAnnoyMnFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdAnnoyMnFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdAnnoyMnFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdAnnoyMnPermImp$AbsVars <- resultsOutAbs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
  coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsVarsConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsVarsConPermimp.pdf")
}

# Plot values disregarding ambient LAeq

pBar <- ggplot(slice(resultsOutAbs.conimp, 1:n()-1)) + geom_col(aes(x=factor(rownames(slice(resultsOutAbs.conimp, 1:n()-1)), levels=rownames(slice(resultsOutAbs.conimp, 1:n()-1))), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimp.pdf")
}

# Plot only positive values disregarding ambient LAeq

resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(slice(resultsOutAbs.conimpPtv, 1:n()-1)) + geom_col(aes(x=factor(rownames(slice(resultsOutAbs.conimpPtv, 1:n()-1)), levels=rownames(slice(resultsOutAbs.conimpPtv, 1:n()-1))), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.svg")
  
  ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.pdf")
}

Selected metric

absVar <- "UASLoudECMAPowAvgBin"

SQM analysis

Replace ambient only stimuli

Here, the ‘ambient only’ stimuli are replaced, as the SQM analysis can incorporate 0 values (whereas dB metrics cannot be handled in this way).

stimDataANum <- rbind(stimDataANum, stimDataA[stimDataA['Stimulus']
                                              == "A1",
                                              colnames(stimDataANum)],
                      stimDataA[stimDataA['Stimulus']
                                              == "A2",
                                              colnames(stimDataANum)])
stimDataANum <- arrange(stimDataANum, Stimulus)

Individual SQMs

Sharpness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASSharpAurISO3PowAvgMaxLR", "UASSharpAurISO305ExMaxLR", "UASSharpAurSHMPowAvgMaxLR", "UASSharpAurSHM05ExMaxLR", "UASSharpAurISO1PowAvgMaxLR", "UASSharpAurISO105ExMaxLR", "UASSharpvBISO1PowAvgMaxLR", "UASSharpvBISO105ExMaxLR", "UASSharpDINPowAvgMaxLR", "UASSharpDIN05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(7041, 905, 4984651, 6513213, 120651)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.6033145

resultsOutSharp$OOB_MAE

[1] 0.4893915

resultsOutSharp$Rsquared

[1] 0.8146039

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.5917107

resultsOutSharp$OOB_MAE

[1] 0.4752747

resultsOutSharp$Rsquared

[1] 0.8242231

# store results
resdAnnoyMnFit['Abs sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdAnnoyMnFit['Abs sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdAnnoyMnFit['Abs sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdAnnoyMnPermImp$AbsSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Sharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnSharpConPermimp.svg", width=8, height=4.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnSharpConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnSharpConPermimp.pdf", width=8, height=4.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnSharpConPermimp.pdf")
}

Selected metric

sharpVar <- "UASSharpAurISO305ExMaxLR"

Tonal loudness and tonality

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(540, 104798, 456464, 87331, 94564)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.6390897

resultsOutTonal1$OOB_MAE

[1] 0.5098136

resultsOutTonal1$Rsquared

[1] 0.7782498

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.6330849

resultsOutTonal1$OOB_MAE

[1] 0.5045541

resultsOutTonal1$Rsquared

[1] 0.7848499

# store results
resdAnnoyMnFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnTonalLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "UASTonLdECMAPowAvgBin"

Tonality without tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(156089, 5860, 10528, 89541, 4685146)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.6030948

resultsOutTonal2$OOB_MAE

[1] 0.4825113

resultsOutTonal2$Rsquared

[1] 0.8186149

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.6032431

resultsOutTonal2$OOB_MAE

[1] 0.4815779

resultsOutTonal2$Rsquared

[1] 0.8154112

# store results
resdAnnoyMnFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnTonalConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "UASTonalSHMInt05ExMaxLR"

Fluctuation strength

Set variables

# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(25107, 546098, 195, 5937, 102658)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.5843817

resultsOutFluct$OOB_MAE

[1] 0.4552687

resultsOutFluct$Rsquared

[1] 0.8202794

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.5847776

resultsOutFluct$OOB_MAE

[1] 0.456084

resultsOutFluct$Rsquared

[1] 0.8194579

# store results
resdAnnoyMnFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImp$AbsFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnFluctConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "UASFluctOV10ExMaxLR"

Roughness

Set variables

# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(4701, 52187, 16589, 65217, 16893)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.5978974

resultsOutRough$OOB_MAE

[1] 0.465474

resultsOutRough$Rsquared

[1] 0.8040179

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.5997308

resultsOutRough$OOB_MAE

[1] 0.4671883

resultsOutRough$Rsquared

[1] 0.8028445

# store results
resdAnnoyMnFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImp$AbsRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnRoughConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnRoughConPermimp.pdf")
}

Selected metric

roughVar <- "UASRoughFZ05ExMaxLR"

Impulsiveness

Set variables

# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(8495, 59867, 5416, 9843, 86)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.6105395

resultsOutImpuls$OOB_MAE

[1] 0.4733032

resultsOutImpuls$Rsquared

[1] 0.7863911

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.603079

resultsOutImpuls$OOB_MAE

[1] 0.4665693

resultsOutImpuls$Rsquared

[1] 0.7960339

# store results
resdAnnoyMnFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnImpulsConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "UASImpulsSHMAvgMaxLR"

SQM and loudness comparison

Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.

Include tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"

seeds <- c(70498, 4, 14986, 453, 864)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.6099365

resultsOutSQMs1$OOB_MAE

[1] 0.4869588

resultsOutSQMs1$Rsquared

[1] 0.8025851

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.60767

resultsOutSQMs1$OOB_MAE

[1] 0.4855243

resultsOutSQMs1$Rsquared

[1] 0.805056

# store results
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"

seeds <- c(546, 57203, 270835, 60592, 8094)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.5870765

resultsOutSQMs2$OOB_MAE

[1] 0.4689816

resultsOutSQMs2$Rsquared

[1] 0.8254522

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.5880773

resultsOutSQMs2$OOB_MAE

[1] 0.469515

resultsOutSQMs2$Rsquared

[1] 0.8255209

# store results
resdAnnoyMnFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf")
}

Difference variables

Next, the difference metrics are analysed,

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)

Set variables

iVars <- c(names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dAnnoyMean"

seeds <- c(568392, 498, 4089, 78132, 741809)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 0.4308342

resultsOutDiffs$OOB_MAE

[1] 0.3567504

resultsOutDiffs$Rsquared

[1] 0.8740203

Train multiple seeds model

resultsOutDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 0.4314902

resultsOutDiffs$OOB_MAE

[1] 0.3566699

resultsOutDiffs$Rsquared

[1] 0.8741682

# store results
resdAnnoyMnFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdAnnoyMnFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdAnnoyMnFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdAnnoyMnPermImp$DiffVars <- resultsOutDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnDiffVarsConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnDiffVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > 0)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnDiffVarsConPermimpPtv.svg")
  
  ggsave(filename="PtAdAnnoyMnDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnDiffVarsConPermimpPtv.pdf")
}

Selected metric

diffVar <- "EPNLLAF50diff"

dSQM analysis

The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.

Individual SQMs

dSharpness

Set variables

iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(84194, 905, 64815, 928054, 625091)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.4523884

resultsOutSharp$OOB_MAE

[1] 0.3595752

resultsOutSharp$Rsquared

[1] 0.8627182

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 0.4536154

resultsOutSharp$OOB_MAE

[1] 0.3606428

resultsOutSharp$Rsquared

[1] 0.8618904

# store results
resdAnnoyMnFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdAnnoyMnFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdAnnoyMnFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdAnnoyMnPermImp$DiffSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMndSharpConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMndSharpConPermimp.pdf")
}

Selected metric

dSharpVar <- "dSharpAurISO3PowAvgMaxLR"

dTonal loudness and dtonality

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "dAnnoyMean"

seeds <- c(561684, 104798, 1536, 48, 48561)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.4541254

resultsOutTonal1$OOB_MAE

[1] 0.3676218

resultsOutTonal1$Rsquared

[1] 0.8616272

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 0.4527969

resultsOutTonal1$OOB_MAE

[1] 0.3685465

resultsOutTonal1$Rsquared

[1] 0.8625633

# store results
resdAnnoyMnFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.6))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMndTonalLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMndTonalLdConPermimp.pdf")
}

Selected metric

dTonLdVar <- "dTonLdECMAPowAvgBin"

dTonality without dtonal loudness

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(410865, 2954, 70812, 203, 7984)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.462319

resultsOutTonal2$OOB_MAE

[1] 0.3667896

resultsOutTonal2$Rsquared

[1] 0.8545831

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 0.4619735

resultsOutTonal2$OOB_MAE

[1] 0.3664541

resultsOutTonal2$Rsquared

[1] 0.8548939

# store results
resdAnnoyMnFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.6))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMndTonalConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMndTonalConPermimp.pdf")
}

Selectec metric

dTonalVar <- "dTonalSHMIntAvgMaxLR"

dFluctuation strength

Set variables

# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(418657, 84, 1630, 18659, 3687)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.4693879

resultsOutFluct$OOB_MAE

[1] 0.3674453

resultsOutFluct$Rsquared

[1] 0.8512241

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 0.4696581

resultsOutFluct$OOB_MAE

[1] 0.368869

resultsOutFluct$Rsquared

[1] 0.8511147

# store results
resdAnnoyMnFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImp$DiffFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnFluctConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnFluctConPermimp.pdf")
}

Selected metric

dFluctVar <- "dFluctOV10ExMaxLR"

dRoughness

Set variables

# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(69851, 85109, 410986, 1563, 896)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.4931997

resultsOutRough$OOB_MAE

[1] 0.372368

resultsOutRough$Rsquared

[1] 0.839582

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 0.4955747

resultsOutRough$OOB_MAE

[1] 0.3727289

resultsOutRough$Rsquared

[1] 0.8381762

# store results
resdAnnoyMnFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImp$DiffRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMndRoughConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMndRoughConPermimp.pdf")
}

Selected metric

dRoughVar <- "dRoughECMA05ExBin"

dImpulsiveness

Set variables

# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "dAnnoyMean"

seeds <- c(418659, 7805, 38475, 65834, 1653)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.5242683

resultsOutImpuls$OOB_MAE

[1] 0.4079044

resultsOutImpuls$Rsquared

[1] 0.8135546

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 0.5207893

resultsOutImpuls$OOB_MAE

[1] 0.4047467

resultsOutImpuls$Rsquared

[1] 0.8173652

# store results
resdAnnoyMnFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMndImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMndImpulsConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMndImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMndImpulsConPermimp.pdf")
}

Selected metric

dImpulsVar <- "dImpulsSHM05ExMaxLR"

dSQM and loudness comparison

Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.

Include dtonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"

seeds <- c(98465, 54163, 6541, 36485, 849675)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.4579162

resultsOutSQMs1$OOB_MAE

[1] 0.3597023

resultsOutSQMs1$Rsquared

[1] 0.8618785

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 0.460312

resultsOutSQMs1$OOB_MAE

[1] 0.363146

resultsOutSQMs1$Rsquared

[1] 0.8597486

# store results
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnDiffSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnDiffSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"

seeds <- c(49865, 7852, 845961, 410583, 36748)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.458182

resultsOutSQMs2$OOB_MAE

[1] 0.360887

resultsOutSQMs2$Rsquared

[1] 0.862422

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 0.4569997

resultsOutSQMs2$OOB_MAE

[1] 0.3596057

resultsOutSQMs2$Rsquared

[1] 0.863165

# store results
resdAnnoyMnFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar

if (saveplots){
  ggsave(filename="PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf")
}

Save the results outputs to file

if (savedata){
  utils::write.csv(resdAnnoyMnFit, paste(outDataPath, "\\ptACRFdAnnoyMnOOBFit.csv", sep=""))
  ii <- 0
  temp = list()
  for (res in resdAnnoyMnPermImp){
    ii <- ii + 1
    temp[[ii]] <- as.data.frame(resdAnnoyMnPermImp[ii])
    names(temp[[ii]]) <- names(resdAnnoyMnPermImp[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFdAnnoyMnConPermimp.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}

(Change to) High annoyance

Initialise results output variables

resdHiAnnoyFit <- data.frame(RMSE = numeric(),
                             MAE = numeric(),
                             Rsquared = numeric())
resdHiAnnoyPermImp <- list()

All variables (absolute and difference)

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)

Set variables

iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
           names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dHighAnnoyPc"

seeds <- c(2, 312, 1897, 465978, 821659)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
              ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutAbsDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 5.468804

resultsOutAbsDiffs$OOB_MAE

[1] 4.327614

resultsOutAbsDiffs$Rsquared

[1] 0.6347201

Train multiple seeds model

resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 5.504329

resultsOutAbsDiffs$OOB_MAE

[1] 4.353992

resultsOutAbsDiffs$Rsquared

[1] 0.6301227

# store results
resdHiAnnoyFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdHiAnnoyFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdHiAnnoyFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdHiAnnoyPermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAllVarsConPermimp.svg")

  ggsave(filename="PtAdHiAnnoyAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAllVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAllVarsConPermimpPtv.svg", width=8, height=12, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAllVarsConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyAllVarsConPermimpPtv.pdf", width=8, height=12, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAllVarsConPermimp.pdf")
}

Absolute variables

Set variables

iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "dHighAnnoyPc"

seeds <- c(578312, 544, 84894, 54654, 153157)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAbsVarsHyperTune.svg")

  ggsave(filename="PtAdHiAnnoyAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAbsVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 20

resultsOutAbs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 5.451148

resultsOutAbs$OOB_MAE

[1] 4.24412

resultsOutAbs$Rsquared

[1] 0.6356777

Train multiple seeds model

resultsOutAbs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 5.450724

resultsOutAbs$OOB_MAE

[1] 4.258414

resultsOutAbs$Rsquared

[1] 0.6356609

# store results
resdHiAnnoyFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdHiAnnoyFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdHiAnnoyFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdHiAnnoyPermImp$AbsVars <- resultsOutAbs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
  coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAbsVarsConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAbsVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimpPtv,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAbsVarsConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAbsVarsConPermimpPtv.svg")
  
  ggsave(filename="PtAdHiAnnoyAbsVarsConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAbsVarsConPermimpPtv.pdf")
}

Selected metric

absVar <- "UASLoudISO3PowAvgBin"

SQM analysis

Replace ambient only stimuli

Here, the ‘ambient only’ stimuli are replaced, as the SQM analysis can incorporate 0 values (whereas dB metrics cannot be handled in this way).

stimDataANum <- rbind(stimDataANum, stimDataA[stimDataA['Stimulus']
                                              == "A1",
                                              colnames(stimDataANum)],
                      stimDataA[stimDataA['Stimulus']
                                              == "A2",
                                              colnames(stimDataANum)])
stimDataANum <- arrange(stimDataANum, Stimulus)

Individual SQMs

Sharpness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASSharpAurISO3PowAvgMaxLR", "UASSharpAurISO305ExMaxLR", "UASSharpAurSHMPowAvgMaxLR", "UASSharpAurSHM05ExMaxLR", "UASSharpAurISO1PowAvgMaxLR", "UASSharpAurISO105ExMaxLR", "UASSharpvBISO1PowAvgMaxLR", "UASSharpvBISO105ExMaxLR", "UASSharpDINPowAvgMaxLR", "UASSharpDIN05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(7041, 905, 4984651, 6513213, 120651)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/2.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.232858

resultsOutSharp$OOB_MAE

[1] 4.190305

resultsOutSharp$Rsquared

[1] 0.6681219

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.241474

resultsOutSharp$OOB_MAE

[1] 4.168303

resultsOutSharp$Rsquared

[1] 0.6679714

# store results
resdHiAnnoyFit['Abs sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdHiAnnoyFit['Abs sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdHiAnnoyFit['Abs sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdHiAnnoyPermImp$AbsSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Sharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoySharpConPermimp.svg", width=8, height=4.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoySharpConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoySharpConPermimp.pdf", width=8, height=4.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoySharpConPermimp.pdf")
}

Selected metric

sharpVar <- "UASSharpAurISO305ExMaxLR"

Tonal loudness and tonality

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(540, 104798, 456464, 87331, 94564)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 5.49583

resultsOutTonal1$OOB_MAE

[1] 4.367014

resultsOutTonal1$Rsquared

[1] 0.6298765

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 5.456376

resultsOutTonal1$OOB_MAE

[1] 4.328807

resultsOutTonal1$Rsquared

[1] 0.6351314

# store results
resdHiAnnoyFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 70))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyTonalLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "UASTonLdECMAPowAvgBin"

Tonality without tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(156089, 5860, 10528, 89541, 4685146)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
                           ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm,
                           minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 5.450319

resultsOutTonal2$OOB_MAE

[1] 4.329233

resultsOutTonal2$Rsquared

[1] 0.6372464

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 5.438749

resultsOutTonal2$OOB_MAE

[1] 4.317752

resultsOutTonal2$Rsquared

[1] 0.6388582

# store results
resdHiAnnoyFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 70))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyTonalConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "UASTonalSHMInt05ExMaxLR"

Fluctuation strength

Set variables

# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(25107, 546098, 195, 5937, 102658)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
                          ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres,
                          nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 5.663893

resultsOutFluct$OOB_MAE

[1] 4.453616

resultsOutFluct$Rsquared

[1] 0.6119349

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 5.631703

resultsOutFluct$OOB_MAE

[1] 4.416416

resultsOutFluct$Rsquared

[1] 0.6156464

# store results
resdHiAnnoyFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImp$AbsFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyFluctConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "UASFluctOV10ExMaxLR"

Roughness

Set variables

# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(4701, 52187, 16589, 65217, 16893)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 5.582484

resultsOutRough$OOB_MAE

[1] 4.337246

resultsOutRough$Rsquared

[1] 0.6209775

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 5.583989

resultsOutRough$OOB_MAE

[1] 4.331894

resultsOutRough$Rsquared

[1] 0.621023

# store results
resdHiAnnoyFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImp$AbsRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyRoughConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyRoughConPermimp.pdf")
}

Selected metric

roughVar <- "UASRoughFZ05ExMaxLR"

Impulsiveness

Set variables

# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(8495, 59867, 5416, 9843, 86)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 5.494835

resultsOutImpuls$OOB_MAE

[1] 4.311435

resultsOutImpuls$Rsquared

[1] 0.6334893

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 5.482274

resultsOutImpuls$OOB_MAE

[1] 4.301147

resultsOutImpuls$Rsquared

[1] 0.6348031

# store results
resdHiAnnoyFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyImpulsConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "UASImpulsSHMAvgMaxLR"

SQM and loudness comparison

Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.

Include tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(70498, 4, 14986, 453, 864)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 5.422008

resultsOutSQMs1$OOB_MAE

[1] 4.293746

resultsOutSQMs1$Rsquared

[1] 0.6411011

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 5.424953

resultsOutSQMs1$OOB_MAE

[1] 4.297558

resultsOutSQMs1$Rsquared

[1] 0.6407472

# store results
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 50))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAbsSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAbsSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(546, 57203, 270835, 60592, 8094)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 5.541448

resultsOutSQMs2$OOB_MAE

[1] 4.390907

resultsOutSQMs2$Rsquared

[1] 0.6273827

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 5.532929

resultsOutSQMs2$OOB_MAE

[1] 4.391987

resultsOutSQMs2$Rsquared

[1] 0.6285529

# store results
resdHiAnnoyFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 50))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf")
}

Difference variables

Next, the difference metrics are analysed,

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)

Set variables

iVars <- c(names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dHighAnnoyPc"

seeds <- c(568392, 498, 4089, 78132, 741809)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/2.25)

Run model

Train preliminary model

nperm <- 30

resultsOutDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 5.719733

resultsOutDiffs$OOB_MAE

[1] 4.428596

resultsOutDiffs$Rsquared

[1] 0.6091806

Train multiple seeds model

resultsOutDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 5.715133

resultsOutDiffs$OOB_MAE

[1] 4.425906

resultsOutDiffs$Rsquared

[1] 0.609705

# store results
resdHiAnnoyFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdHiAnnoyFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdHiAnnoyFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdHiAnnoyPermImp$DiffVars <- resultsOutDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyDiffVarsConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyDiffVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > 0)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyDiffVarsConPermimpPtv.svg")
  
  ggsave(filename="PtAdHiAnnoyDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyDiffVarsConPermimpPtv.pdf")
}

Selected metric

diffVar <- "EPNLLAeqdiff"

dSQM analysis

The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.

Individual SQMs

dSharpness

Set variables

iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(84194, 905, 64815, 928054, 625091, 582031)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.817107

resultsOutSharp$OOB_MAE

[1] 4.431913

resultsOutSharp$Rsquared

[1] 0.5844272

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 5.794601

resultsOutSharp$OOB_MAE

[1] 4.410331

resultsOutSharp$Rsquared

[1] 0.5874878

# store results
resdHiAnnoyFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdHiAnnoyFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdHiAnnoyFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdHiAnnoyPermImp$DiffSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoydSharpConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoydSharpConPermimp.pdf")
}

Selected metric

dSharpVar <- "dSharpAurSHMPowAvgMaxLR"

dTonal loudness and dtonality

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "dHighAnnoyPc"

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

seeds <- c(561684, 104798, 1536, 48, 48561)

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 6.71205

resultsOutTonal1$OOB_MAE

[1] 5.189994

resultsOutTonal1$Rsquared

[1] 0.4460736

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 6.722309

resultsOutTonal1$OOB_MAE

[1] 5.199688

resultsOutTonal1$Rsquared

[1] 0.4443759

# store results
resdHiAnnoyFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 70))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoydTonalLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoydTonalLdConPermimp.pdf")
}

Selected metric

dTonLdVar <- "dTonLdECMAPowAvgBin"

dTonality without dtonal loudness

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(410865, 2954, 70812, 203, 7984)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

NA
NA

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 6.581687

resultsOutTonal2$OOB_MAE

[1] 5.067818

resultsOutTonal2$Rsquared

[1] 0.4671435

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 6.580352

resultsOutTonal2$OOB_MAE

[1] 5.068529

resultsOutTonal2$Rsquared

[1] 0.4673965

# store results
resdHiAnnoyFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 70))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoydTonalConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoydTonalConPermimp.pdf")
}

Selected metric

dTonalVar <- "dTonalSHMIntAvgMaxLR"

dFluctuation strength

Set variables

# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(418657, 84, 1630, 18659, 3687)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 6.714071

resultsOutFluct$OOB_MAE

[1] 5.1887

resultsOutFluct$Rsquared

[1] 0.4455587

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 6.702409

resultsOutFluct$OOB_MAE

[1] 5.191661

resultsOutFluct$Rsquared

[1] 0.4472641

# store results
resdHiAnnoyFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImp$DiffFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyFluctConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyFluctConPermimp.pdf")
}

Selected metric

dFluctVar <- "dFluctSHM10ExBin"

dRoughness

Set variables

# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(69851, 85109, 410986, 1563, 896)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

Warning: Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')

# print model prediction results
resultsOutRough$OOB_RMSE

NULL

resultsOutRough$OOB_MAE

NULL

resultsOutRough$Rsquared

NULL

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 6.950843

resultsOutRough$OOB_MAE

[1] 5.259601

resultsOutRough$Rsquared

[1] 0.4054285

# store results
resdHiAnnoyFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImp$DiffRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoydRoughConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoydRoughConPermimp.pdf")
}

Selected metric

dRoughVar <- "dRoughFZ05ExMaxLR"

dImpulsiveness

Set variables

# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "dHighAnnoyPc"

seeds <- c(418659, 7805, 38475, 65834, 1653)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 6.674278

resultsOutImpuls$OOB_MAE

[1] 5.175702

resultsOutImpuls$Rsquared

[1] 0.4517669

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 6.673098

resultsOutImpuls$OOB_MAE

[1] 5.185144

resultsOutImpuls$Rsquared

[1] 0.4520846

# store results
resdHiAnnoyFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoydImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoydImpulsConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoydImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoydImpulsConPermimp.pdf")
}

Selected metric

dImpulsVar <- "dImpulsSHMAvgMaxLR"

dSQM and loudness comparison

Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.

Include dtonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(98465, 54163, 6541, 36485, 849675)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

Warning: Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')

# print model prediction results
resultsOutSQMs1$OOB_RMSE

NULL

resultsOutSQMs1$OOB_MAE

NULL

resultsOutSQMs1$Rsquared

NULL

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 5.939049

resultsOutSQMs1$OOB_MAE

[1] 4.654886

resultsOutSQMs1$Rsquared

[1] 0.5768752

# store results
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyDiffSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyDiffSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"

seeds <- c(49865, 7852, 845961, 410583, 36748)

Hyperparameter tuning

p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 5.955857

resultsOutSQMs2$OOB_MAE

[1] 4.656295

resultsOutSQMs2$Rsquared

[1] 0.5742314

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 5.966773

resultsOutSQMs2$OOB_MAE

[1] 4.669058

resultsOutSQMs2$Rsquared

[1] 0.5722321

# store results
resdHiAnnoyFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar

if (saveplots){
  ggsave(filename="PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf")
}

Save the results outputs to file

if (savedata){
  utils::write.csv(resdHiAnnoyFit, paste(outDataPath, "\\ptACRFdHiAnnoyOOBFit.csv", sep=""))
  ii <- 0
  temp = list()
  for (res in resdHiAnnoyPermImp){
    ii <- ii + 1
    temp[[ii]] <- as.data.frame(resdHiAnnoyPermImp[ii])
    names(temp[[ii]]) <- names(resdHiAnnoyPermImp[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFdHiAnnoyConPermimp.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}

Noticeability

Initialise results output variables

resNoticeFit <- data.frame(RMSE = numeric(),
                             MAE = numeric(),
                             Rsquared = numeric())
resNoticePermImp <- list()

All variables (absolute and difference)

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimNoticeDataANum <- stimNoticeDataANum[complete.cases(stimNoticeDataANum),]
stimNoticeDataANum$UASLAeq <- as.numeric(stimNoticeDataANum$UASLAeq)
stimNoticeDataANum$SNRlevel <- as.numeric(stimNoticeDataANum$SNRlevel)

Set variables

iVars <- names(stimNoticeDataANum)[which(names(stimNoticeDataANum) == 'UASLAeq'):which(names(stimNoticeDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
           names(stimNoticeDataANum)[which(colnames(stimNoticeDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimNoticeDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "NoticedPcFilt"

seeds <- c(2, 312, 1897, 465978, 821659)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/2.75)

Run model

Train preliminary model

nperm <- 10

resultsOutAbsDiffs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 10.23042

resultsOutAbsDiffs$OOB_MAE

[1] 7.618346

resultsOutAbsDiffs$Rsquared

[1] 0.9097534

Train multiple seeds model

resultsOutAbsDiffs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbsDiffs$OOB_RMSE

[1] 10.19738

resultsOutAbsDiffs$OOB_MAE

[1] 7.652439

resultsOutAbsDiffs$Rsquared

[1] 0.9103752

# store results
resNoticeFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resNoticeFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resNoticeFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resNoticePermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAllVarsConPermimp.svg")

  ggsave(filename="PtANoticeAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAllVarsConPermimp.pdf")
}

# Plot only positive values

resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeAllVarsConPermimpPtv.svg", width=8, height=12, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAllVarsConPermimp.svg")
  
  ggsave(filename="PtANoticeAllVarsConPermimpPtv.pdf", width=8, height=12, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAllVarsConPermimp.pdf")
}

Absolute variables

Set variables

iVars <- names(stimNoticeDataANum)[which(names(stimNoticeDataANum) == 'UASLAeq'):which(names(stimNoticeDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "NoticedPcFilt"

seeds <- c(578312, 544, 84894, 54654, 153157)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

if (saveplots){
  ggsave(filename="PtANoticeAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAbsVarsHyperTune.svg")

  ggsave(filename="PtANoticeAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAbsVarsHyperTune.pdf")
}

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 20

resultsOutAbs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 16.59018

resultsOutAbs$OOB_MAE

[1] 12.64339

resultsOutAbs$Rsquared

[1] 0.7829788

Train multiple seeds model

resultsOutAbs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutAbs$OOB_RMSE

[1] 16.50718

resultsOutAbs$OOB_MAE

[1] 12.55933

resultsOutAbs$Rsquared

[1] 0.7846705

# store results
resNoticeFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resNoticeFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resNoticeFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resNoticePermImp$AbsVars <- resultsOutAbs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
  coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAbsVarsConPermimp.svg")
  
  ggsave(filename="PtANoticeAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAbsVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
                                          rownames_to_column('Metric') |>
                                                filter_if(is.numeric, all_vars(. > 0)) |>
                                                      column_to_rownames('Metric')

pBar <- ggplot(resultsOutAbs.conimpPtv,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeAbsVarsConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAbsVarsConPermimpPtv.svg")
  
  ggsave(filename="PtANoticeAbsVarsConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAbsVarsConPermimpPtv.pdf")
}

Selected metric

absVar <- "UASPNLmaxMaxLR"

SQM analysis

Replace ambient only stimuli

Here, the ‘ambient only’ stimuli are not replaced, as the analysis cannot handle missing values for dB metrics (PNL).

# 
# stimNoticeDataANum <- rbind(stimNoticeDataANum, stimNoticeDataA[stimNoticeDataA['Stimulus']
#                                               == "A1",
#                                               colnames(stimNoticeDataANum)],
#                       stimNoticeDataA[stimNoticeDataA['Stimulus']
#                                               == "A2",
#                                               colnames(stimNoticeDataANum)])
# stimNoticeDataANum <- arrange(stimNoticeDataANum, Stimulus)

Individual SQMs

Sharpness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASSharpAurISO3PowAvgMaxLR", "UASSharpAurISO305ExMaxLR", "UASSharpAurSHMPowAvgMaxLR", "UASSharpAurSHM05ExMaxLR", "UASSharpAurISO1PowAvgMaxLR", "UASSharpAurISO105ExMaxLR", "UASSharpvBISO1PowAvgMaxLR", "UASSharpvBISO105ExMaxLR", "UASSharpDINPowAvgMaxLR", "UASSharpDIN05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(7041, 905, 4984651, 6513213, 120651)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
  
# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 18.90396

resultsOutSharp$OOB_MAE

[1] 14.44931

resultsOutSharp$Rsquared

[1] 0.7027075

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 18.96332

resultsOutSharp$OOB_MAE

[1] 14.46537

resultsOutSharp$Rsquared

[1] 0.6996689

# store results
resNoticeFit['Abs sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resNoticeFit['Abs sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resNoticeFit['Abs sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resNoticePermImp$AbsSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Sharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeSharpConPermimp.svg", width=8, height=4.9, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeSharpConPermimp.svg")
  
  ggsave(filename="PtANoticeSharpConPermimp.pdf", width=8, height=4.9, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeSharpConPermimp.pdf")
}

Selected metric

sharpVar <- "UASSharpAurISO305ExMaxLR"

Tonal loudness and tonality

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(540, 104798, 456464, 87331, 94564)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 18.47209

resultsOutTonal1$OOB_MAE

[1] 14.20136

resultsOutTonal1$Rsquared

[1] 0.7179216

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 18.19492

resultsOutTonal1$OOB_MAE

[1] 14.0021

resultsOutTonal1$Rsquared

[1] 0.7280624

# store results
resNoticeFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resNoticeFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resNoticeFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resNoticePermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 1000))
pBar

if (saveplots){
  ggsave(filename="PtANoticeTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeTonalLdConPermimp.svg")
  
  ggsave(filename="PtANoticeTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "UASTonLdECMA05ExBin"

Tonality without tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(156089, 5860, 10528, 89541, 4685146)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 17.0126

resultsOutTonal2$OOB_MAE

[1] 13.20926

resultsOutTonal2$Rsquared

[1] 0.7729116

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 17.06205

resultsOutTonal2$OOB_MAE

[1] 13.236

resultsOutTonal2$Rsquared

[1] 0.7718614

# store results
resNoticeFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resNoticeFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resNoticeFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resNoticePermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 1000))
pBar

if (saveplots){
  ggsave(filename="PtANoticeTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeTonalConPermimp.svg")
  
  ggsave(filename="PtANoticeTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "UASTonalSHMInt05ExMaxLR"

Fluctuation strength

Set variables

# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(25107, 546098, 195, 5937, 102658)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.5)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 16.91354

resultsOutFluct$OOB_MAE

[1] 12.93498

resultsOutFluct$Rsquared

[1] 0.7866218

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 16.96583

resultsOutFluct$OOB_MAE

[1] 12.82457

resultsOutFluct$Rsquared

[1] 0.7842884

# store results
resNoticeFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resNoticeFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resNoticeFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resNoticePermImp$AbsFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeFluctConPermimp.svg")
  
  ggsave(filename="PtANoticeFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "UASFluctOV10ExMaxLR"

Roughness

Set variables

# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(4701, 52187, 16589, 65217, 16893)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 17.26411

resultsOutRough$OOB_MAE

[1] 13.17981

resultsOutRough$Rsquared

[1] 0.7714465

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 16.75907

resultsOutRough$OOB_MAE

[1] 12.82419

resultsOutRough$Rsquared

[1] 0.7847338

# store results
resNoticeFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resNoticeFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resNoticeFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resNoticePermImp$AbsRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeRoughConPermimp.svg")
  
  ggsave(filename="PtANoticeRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeRoughConPermimp.pdf")
}

Selected metric

roughVar <- "UASRoughFZ05ExMaxLR"

Impulsiveness

Set variables

# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(8495, 59867, 5416, 9843, 86)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 16.24906

resultsOutImpuls$OOB_MAE

[1] 12.47619

resultsOutImpuls$Rsquared

[1] 0.7966846

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 16.26418

resultsOutImpuls$OOB_MAE

[1] 12.46156

resultsOutImpuls$Rsquared

[1] 0.7959352

# store results
resNoticeFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resNoticeFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resNoticeFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resNoticePermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeImpulsConPermimp.svg")
  
  ggsave(filename="PtANoticeImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "UASImpulsSHMAvgMaxLR"

SQM and loudness comparison

Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.

Include tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "NoticedPcFilt"

seeds <- c(70498, 4, 14986, 453, 864)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 18.88386

resultsOutSQMs1$OOB_MAE

[1] 14.50936

resultsOutSQMs1$Rsquared

[1] 0.7014957

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 18.76572

resultsOutSQMs1$OOB_MAE

[1] 14.41058

resultsOutSQMs1$Rsquared

[1] 0.7076868

# store results
resNoticeFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resNoticeFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resNoticeFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resNoticePermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 800))
pBar

if (saveplots){
  ggsave(filename="PtANoticeAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAbsSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtANoticeAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAbsSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "NoticedPcFilt"

seeds <- c(546, 57203, 270835, 60592, 8094)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 18.88129

resultsOutSQMs2$OOB_MAE

[1] 14.46072

resultsOutSQMs2$Rsquared

[1] 0.7037345

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 18.78148

resultsOutSQMs2$OOB_MAE

[1] 14.40389

resultsOutSQMs2$Rsquared

[1] 0.7078282

# store results
resNoticeFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resNoticeFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resNoticeFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resNoticePermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 800))
pBar

if (saveplots){
  ggsave(filename="PtANoticeAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeAbsSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtANoticeAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeAbsSQMsNoTonLdConPermimp.pdf")
}

Difference variables

Next, the difference metrics are analysed,

Remove ambient only stimuli

Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.

stimNoticeDataANum <- stimNoticeDataANum[complete.cases(stimNoticeDataANum),]
stimNoticeDataANum$UASLAeq <- as.numeric(stimNoticeDataANum$UASLAeq)
stimNoticeDataANum$SNRlevel <- as.numeric(stimNoticeDataANum$SNRlevel)

Set variables

iVars <- c(names(stimNoticeDataANum)[which(colnames(stimNoticeDataANum)=="LAeqLAF90diff"):
                               which(colnames(stimNoticeDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "NoticedPcFilt"

seeds <- c(568392, 498, 4089, 78132, 741809)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/3.5)

Run model

Train preliminary model

nperm <- 30

resultsOutDiffs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 10.09164

resultsOutDiffs$OOB_MAE

[1] 7.746637

resultsOutDiffs$Rsquared

[1] 0.9128949

Train multiple seeds model

resultsOutDiffs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutDiffs$OOB_RMSE

[1] 10.15461

resultsOutDiffs$OOB_MAE

[1] 7.781117

resultsOutDiffs$Rsquared

[1] 0.9118708

# store results
resNoticeFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resNoticeFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resNoticeFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resNoticePermImp$DiffVars <- resultsOutDiffs$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeDiffVarsConPermimp.svg")
  
  ggsave(filename="PtANoticeDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeDiffVarsConPermimp.pdf")
}

# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
                                        rownames_to_column('Metric') |>
                                            filter_if(is.numeric, all_vars(. > 0)) |>
                                                  column_to_rownames('Metric')

pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticeDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeDiffVarsConPermimpPtv.svg")
  
  ggsave(filename="PtANoticeDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeDiffVarsConPermimpPtv.pdf")
}

Selected metric

diffVar <- "Detect0p5dBMaxMaxLR"

dSQM analysis

The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.

Individual SQMs

dSharpness

Set variables

iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(84194, 905, 928054, 625091, 582031)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutSharp <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 10.40193

resultsOutSharp$OOB_MAE

[1] 8.108278

resultsOutSharp$Rsquared

[1] 0.9042801

Train multiple seeds model

resultsOutSharp <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSharp$OOB_RMSE

[1] 10.43319

resultsOutSharp$OOB_MAE

[1] 8.137405

resultsOutSharp$Rsquared

[1] 0.9035948

# store results
resNoticeFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resNoticeFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resNoticeFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resNoticePermImp$DiffSharp <- resultsOutSharp$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticedSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtANoticedSharpConPermimp.svg")
  
  ggsave(filename="PtANoticedSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticedSharpConPermimp.pdf")
}

Selected metric

sharpVar <- "dSharpAurISO305ExMaxLR"

dTonal loudness and dtonality

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "NoticedPcFilt"

seeds <- c(561684, 104798, 1536, 48, 48561)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <-1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality with tonal loudness

nperm <- 5

resultsOutTonal1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 10.46408

resultsOutTonal1$OOB_MAE

[1] 8.013044

resultsOutTonal1$Rsquared

[1] 0.9021881

Train multiple seeds model

# Tonality with tonal loudness

resultsOutTonal1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal1$OOB_RMSE

[1] 10.45256

resultsOutTonal1$OOB_MAE

[1] 8.001321

resultsOutTonal1$Rsquared

[1] 0.9026108

# store results
resNoticeFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resNoticeFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resNoticeFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resNoticePermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1200))
pBar

if (saveplots){
  ggsave(filename="PtANoticedTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
  unlink("PtANoticedTonalLdConPermimp.svg")
  
  ggsave(filename="PtANoticedTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticedTonalLdConPermimp.pdf")
}

Selected metric

tonLdVar <- "dTonLdECMAPowAvgBin"

dTonality without dtonal loudness

Set variables

iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(410865, 2954, 70812, 203, 7984)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

# Tonality

nperm <- 5

resultsOutTonal2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 12.53815

resultsOutTonal2$OOB_MAE

[1] 9.798852

resultsOutTonal2$Rsquared

[1] 0.8608015

Train multiple seeds model

# Tonality

resultsOutTonal2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutTonal2$OOB_RMSE

[1] 12.57501

resultsOutTonal2$OOB_MAE

[1] 9.807799

resultsOutTonal2$Rsquared

[1] 0.8595653

# store results
resNoticeFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resNoticeFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resNoticeFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resNoticePermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1200))
pBar

if (saveplots){
  ggsave(filename="PtANoticedTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticedTonalConPermimp.svg")
  
  ggsave(filename="PtANoticedTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticedTonalConPermimp.pdf")
}

Selected metric

tonalVar <- "dTonalECMAAvgMaxLR"

dFluctuation strength

Set variables

# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(418657, 84, 1630, 18659, 3687)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutFluct <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 11.38464

resultsOutFluct$OOB_MAE

[1] 8.885649

resultsOutFluct$Rsquared

[1] 0.8846925

Train multiple seeds model

resultsOutFluct <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutFluct$OOB_RMSE

[1] 11.39714

resultsOutFluct$OOB_MAE

[1] 8.868102

resultsOutFluct$Rsquared

[1] 0.8843449

# store results
resNoticeFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resNoticeFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resNoticeFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resNoticePermImp$DiffFluct <- resultsOutFluct$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticedFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeFluctConPermimp.svg")
  
  ggsave(filename="PtANoticedFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeFluctConPermimp.pdf")
}

Selected metric

fluctVar <- "dFluctSHM05ExBin"

dRoughness

Set variables

# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(69851, 85109, 410986, 1563, 896)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutRough <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 12.37658

resultsOutRough$OOB_MAE

[1] 9.87597

resultsOutRough$Rsquared

[1] 0.8723291

Train multiple seeds model

resultsOutRough <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutRough$OOB_RMSE

[1] 12.37394

resultsOutRough$OOB_MAE

[1] 9.849881

resultsOutRough$Rsquared

[1] 0.8718619

# store results
resNoticeFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resNoticeFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resNoticeFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resNoticePermImp$DiffRough <- resultsOutRough$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticedRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticedRoughConPermimp.svg")
  
  ggsave(filename="PtANoticedRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticedRoughConPermimp.pdf")
}

Selected metric

roughVar <- "dRoughECMA05ExBin"

dImpulsiveness

Set variables

# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "NoticedPcFilt"

seeds <- c(418659, 7805, 38475, 65834, 1653)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 251
mtry <- as.integer(length(iVars)/1.25)

Run model

Train preliminary model

nperm <- 5

resultsOutImpuls <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 11.85234

resultsOutImpuls$OOB_MAE

[1] 9.182759

resultsOutImpuls$Rsquared

[1] 0.8773988

Train multiple seeds model

resultsOutImpuls <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutImpuls$OOB_RMSE

[1] 11.92896

resultsOutImpuls$OOB_MAE

[1] 9.255452

resultsOutImpuls$Rsquared

[1] 0.8757174

# store results
resNoticeFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resNoticeFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resNoticeFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resNoticePermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar

if (saveplots){
  ggsave(filename="PtANoticedImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
  unlink("PtANoticedImpulsConPermimp.svg")
  
  ggsave(filename="PtANoticedImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticedImpulsConPermimp.pdf")
}

Selected metric

impulsVar <- "dImpulsSHM05ExMaxLR"

dSQM and loudness comparison

Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.

Include dtonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "NoticedPcFilt"

seeds <- c(98465, 54163, 6541, 36485, 849675)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 501
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 10.79004

resultsOutSQMs1$OOB_MAE

[1] 8.514653

resultsOutSQMs1$Rsquared

[1] 0.9054054

Train multiple seeds model

resultsOutSQMs1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs1$OOB_RMSE

[1] 10.5503

resultsOutSQMs1$OOB_MAE

[1] 8.348676

resultsOutSQMs1$Rsquared

[1] 0.9084059

# store results
resNoticeFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resNoticeFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resNoticeFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resNoticePermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 600))
pBar

if (saveplots){
  ggsave(filename="PtANoticeDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeDiffSQMsTonLdConPermimp.svg")
  
  ggsave(filename="PtANoticeDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeDiffSQMsTonLdConPermimp.pdf")
}

Exclude tonal loudness

Set variables

iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "NoticedPcFilt"

seeds <- c(49865, 7852, 845961, 410583, 36748)

Hyperparameter tuning

p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
             ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)

p

Selected hyperparameters

ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)

Run model

Train preliminary model

nperm <- 5

resultsOutSQMs2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 11.8702

resultsOutSQMs2$OOB_MAE

[1] 9.243409

resultsOutSQMs2$Rsquared

[1] 0.8859168

Train multiple seeds model

resultsOutSQMs2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)

# print model prediction results
resultsOutSQMs2$OOB_RMSE

[1] 11.81577

resultsOutSQMs2$OOB_MAE

[1] 9.205458

resultsOutSQMs2$Rsquared

[1] 0.8863226

# store results
resNoticeFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resNoticeFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resNoticeFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resNoticePermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp

Plot results

par(mai=c(0,3,0,0))

# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))

pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 600))
pBar

if (saveplots){
  ggsave(filename="PtANoticeDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
  unlink("PtANoticeDiffSQMsNoTonLdConPermimp.svg")
  
  ggsave(filename="PtANoticeDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
  unlink("PtANoticeDiffSQMsNoTonLdConPermimp.pdf")
}

Save the results outputs to file

if (savedata){
  utils::write.csv(resNoticeFit, paste(outDataPath, "\\ptACRFNoticeOOBFit.csv", sep=""))
  ii <- 0
  temp = list()
  for (res in resNoticePermImp){
    ii <- ii + 1
    temp[[ii]] <- as.data.frame(resNoticePermImp[ii])
    names(temp[[ii]]) <- names(resNoticePermImp[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFNoticeConPermimp.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}

Part A summary

Summary of results for Part A

With tonal loudness

Absolute variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpTbl <- as.data.frame(resdAnnoyMnPermImp$AbsSQMs1/max(resdAnnoyMnPermImp$AbsSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpTbl)[2] <- "dAnnoy"

resdHiAnnoyAbsPermImpTbl <- as.data.frame(resdHiAnnoyPermImp$AbsSQMs1/max(resdHiAnnoyPermImp$AbsSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpTbl)[2] <- "dHiAnnoy"

resNoticeAbsPermImpTbl <- as.data.frame(resNoticePermImp$AbsSQMs1/max(resNoticePermImp$AbsSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resNoticeAbsPermImpTbl)[2] <- "dNotice"

# merge the dataframes
resAbsPermImpTbl <- list(resdAnnoyMnAbsPermImpTbl, resdHiAnnoyAbsPermImpTbl, resNoticeAbsPermImpTbl) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resAbsPermImpTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resAbsPermImpTbl[is.na(resAbsPermImpTbl)] <- 0

resAbs <- tidyr::pivot_longer(resAbsPermImpTbl, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbs <- resAbs |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resAbs$Outcome <- factor(resAbs$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbs) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35',  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resAbs$Variable)))
pBar

if (saveplots){
  ggsave(filename="PtAcrfAbsSQMsSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAcrfAbsSQMsSummary.svg")
  
  ggsave(filename="PtAcrfAbsSQMsSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAcrfAbsSQMsSummary.pdf")
}

Difference variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpTbl <- as.data.frame(resdAnnoyMnPermImp$DiffSQMs1/max(resdAnnoyMnPermImp$DiffSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpTbl)[2] <- "dAnnoy"

resdHiAnnoyDiffPermImpTbl <- as.data.frame(resdHiAnnoyPermImp$DiffSQMs1/max(resdHiAnnoyPermImp$DiffSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpTbl)[2] <- "dHiAnnoy"

resNoticeDiffPermImpTbl <- as.data.frame(resNoticePermImp$DiffSQMs1/max(resNoticePermImp$DiffSQMs1)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resNoticeDiffPermImpTbl)[2] <- "dNotice"

# merge the dataframes
resDiffPermImpTbl <- list(resdAnnoyMnDiffPermImpTbl, resdHiAnnoyDiffPermImpTbl, resNoticeDiffPermImpTbl) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resDiffPermImpTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resDiffPermImpTbl[is.na(resDiffPermImpTbl)] <- 0

resDiff <- tidyr::pivot_longer(resDiffPermImpTbl, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiff <- resDiff |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resDiff$Outcome <- factor(resDiff$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiff) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35',  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiff$Variable)))
pBar

if (saveplots){
  ggsave(filename="PtAcrfDiffSQMsSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAcrfDiffSQMsSummary.svg")
  
  ggsave(filename="PtAcrfDiffSQMsSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAcrfDiffSQMsSummary.pdf")
}

No tonal loudness

Absolute variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpNoTonLdTbl <- as.data.frame(resdAnnoyMnPermImp$AbsSQMs2/max(resdAnnoyMnPermImp$AbsSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpNoTonLdTbl)[2] <- "dAnnoy"

resdHiAnnoyAbsPermImpNoTonLdTbl <- as.data.frame(resdHiAnnoyPermImp$AbsSQMs2/max(resdHiAnnoyPermImp$AbsSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpNoTonLdTbl)[2] <- "dHiAnnoy"

resNoticeAbsPermImpNoTonLdTbl <- as.data.frame(resNoticePermImp$AbsSQMs2/max(resNoticePermImp$AbsSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resNoticeAbsPermImpNoTonLdTbl)[2] <- "dNotice"

# merge the dataframes
resAbsNoTonLdPermImpNoTonLdTbl <- list(resdAnnoyMnAbsPermImpNoTonLdTbl, resdHiAnnoyAbsPermImpNoTonLdTbl, resNoticeAbsPermImpNoTonLdTbl) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resAbsNoTonLdPermImpNoTonLdTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resAbsNoTonLdPermImpNoTonLdTbl[is.na(resAbsNoTonLdPermImpNoTonLdTbl)] <- 0

resAbsNoTonLd <- tidyr::pivot_longer(resAbsNoTonLdPermImpNoTonLdTbl, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbsNoTonLd <- resAbsNoTonLd |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resAbsNoTonLd$Outcome <- factor(resAbsNoTonLd$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbsNoTonLd) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35',  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resAbsNoTonLd$Variable)))
pBar

if (saveplots){
  ggsave(filename="PtAcrfAbsSQMsNoTonLdSummary.svg", width=8, height=3.7, path=file.path(outFigPath, "svg"))
  unlink("PtAcrfAbsSQMsNoTonLdSummary.svg")
  
  ggsave(filename="PtAcrfAbsSQMsNoTonLdSummary.pdf", width=8, height=3.7, path=file.path(outFigPath, "pdf"))
  unlink("PtAcrfAbsSQMsNoTonLdSummary.pdf")
}

Difference variables

# combine the annoyance perm importance results

# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpNoTonLdTbl <- as.data.frame(resdAnnoyMnPermImp$DiffSQMs2/max(resdAnnoyMnPermImp$DiffSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpNoTonLdTbl)[2] <- "dAnnoy"

resdHiAnnoyDiffPermImpNoTonLdTbl <- as.data.frame(resdHiAnnoyPermImp$DiffSQMs2/max(resdHiAnnoyPermImp$DiffSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpNoTonLdTbl)[2] <- "dHiAnnoy"

resNoticeDiffPermImpNoTonLdTbl <- as.data.frame(resNoticePermImp$DiffSQMs2/max(resNoticePermImp$DiffSQMs2)) |>
  tibble::rownames_to_column(var='Variable')
colnames(resNoticeDiffPermImpNoTonLdTbl)[2] <- "dNotice"

# merge the dataframes
resDiffNoTonLdPermImpNoTonLdTbl <- list(resdAnnoyMnDiffPermImpNoTonLdTbl, resdHiAnnoyDiffPermImpNoTonLdTbl, resNoticeDiffPermImpNoTonLdTbl) |>
  purrr::reduce(merge, by = c('Variable'), all = T)

# rename the columns
colnames(resDiffNoTonLdPermImpNoTonLdTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resDiffNoTonLdPermImpNoTonLdTbl[is.na(resDiffNoTonLdPermImpNoTonLdTbl)] <- 0

resDiffNoTonLd <- tidyr::pivot_longer(resDiffNoTonLdPermImpNoTonLdTbl, cols=-Variable, names_to="Outcome", values_to="Imp")

# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiffNoTonLd <- resDiffNoTonLd |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
   mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))

# Reorder outcome levels
resDiffNoTonLd$Outcome <- factor(resDiffNoTonLd$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))

# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiffNoTonLd) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35',  width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiffNoTonLd$Variable)))
pBar

if (saveplots){
  ggsave(filename="PtAcrfDiffSQMsNoTonLdSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
  unlink("PtAcrfDiffSQMsNoTonLdSummary.svg")
  
  ggsave(filename="PtAcrfDiffSQMsNoTonLdSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
  unlink("PtAcrfDiffSQMsNoTonLdSummary.pdf")
}

Save the results outputs to file

# Make a list of the summary results
resSummary <- list(resAbs, resDiff, resAbsNoTonLd, resDiffNoTonLd)

# Save the results
if (savedata){
  ii <- 0
  temp = list()
  for (res in resSummary){
    ii <- ii + 1
    temp[[ii]] <- data.frame(resSummary[ii])
  }
  openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFSummary.xlsx",
                                   sep=""),
                       rowNames=TRUE)
}